General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 



National Research 
Council Canada 


Conseil national 
de recherches Canada 



A NUMERICAL DETERMINATION 
OF THE BOW SHOCK WAVE IN 
TRANSONIC AXISYMMETRIC 
FLOW ABOUT BLUNT BODIES 


ASA-TM-X-72448) A KOMEEICAL DETERMINATION 
THE BOH SHOCK HAVE IN TRANSONIC 
ISYMMETRIC FLOW AEODT BLUNT BODIES (NASA) 
p HC $4.25 CSCL 01/1 


03/02 


N75-29035 


Unclas 

31flfl‘i 


BY 


D.J. JONES. J.C. SOUTH. JR. 
NATIONAL AERONAUTICAL ESTABLISHMENT 


OTTAWA 
MAY 1975 


NRC IMO. 14765 
ISSN 0077 5541 


AERONAUTICAL 
REPORT 
LR 586 


A NUMERICAL DETERMINATION OF THE BOW SHOCK WAVE 
IN TRANSONIC AXISYMMETRIC FLOW ABOUT BLUNT BODIES 


UNE DETERMINATION NUMERIQUE DE L’ONDE DE CHOC DE TETE 
AUTOUR DE CORPS EMOUSSES DANS UN ECOULEMENT TRANSONIQUE 

AXISYMETRIQUE 


by /par 

D.J. JONES*, J.C. SOUTH, JR.f 


* Part of this work was carried out at NASA Langley where this author was on leave November 1973 — 
May 1974, 
t NASA Langley 


L.H. Ohman, Head /Chef 

High Speed Aerodynamics Laboratory/ 

Laboratoire d’aerodynamique a hautes vitesses 


F.R. Thurston 
Director /Directeur 



SUMMARY 


A numerical method is developed for calculating axisymmetric 
transonic (M > 1) flow about a blunt body. The bow shock wave location 
is of particular interest. A Rankine Hugoniot jump is applied at the shock 
while relaxation on the isentropic equation of motion is used between shock 
and body. The shock wave is adjusted by a Newton type iteration scheme. 
Results are given for a sphere in the Mach number range 1.62 down to 1.02. 


RESUME 


Une methode numerique est develop pee pour calculer I’ecoulement 
transonique axisymetrique (M > 1) autour d’un corps emousse. Le lieu et la 
forme de I’onde de tete revet un int^ret particulier. Les equations de Rankine 
Hugoniot sont utilisees a la traversee de Tonde de tete et une methode de 
relaxation est appliquee a I’equation isentropique du mouvement dans la 
couche de choc. La position du choc est reglee par la methode iterative de 
Newton. Les resultats pour une sphere sont presentes pour des nombres de 
Mach allant de 1.62 a 1.02. 


(iii) 
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A NUMERICAL DETERMINATION OF THE BOW SHOCK WAVE 
IN TRANSONIC AXISYMMETRIC FLOW ABOUT BLUNT BODIES 


1.0 INTRODUCTION 

For many years the problem of determining the bovir shock wave location has been tackled 
by various authors. Several successful methods have resulted including the method of integral relations 
(Schemes I, II and III) of Belotserkovskiy (Ref. 4) and time dependent approaches such as Barnwell’s 
(Ref. 3). 


Belotserkovskiy used Schemes I and II, in which one of the independent variables is made 
discrete while the other is kept continuous, to compute solutions for flow about spheres and other 
shapes and was successful in computing flow characteristics for Mach numbers ranging from hypersonic 
down to about 1.5. For Mach numbers lower than this convergence was difficult to obtain so he de- 
vised Scheme III in which both independent variables were discretized and dependent variables were 
represented by polynomials. In this way solutions could be obtained for Mach numbers down to 1.05. 

Barnwell’s method, which in general will treat bodies with discontinuous slope as well a'j 
bodies at incidence, starts from an assumed shock wave and integrates the time dependent equations 
of motion xintil a steady state is reached. Solutions can be obtained efficiently for Mach numbers as 
low as about 1.3 but below this the method may produce kinky shocks or else is very time consuming. 
For example another time dependent method due to Aungier (Ref. 9) took 16 hours to compute the 
flow about a hemisphere cylinder at M = 1.1 on an IBM 370/155 (quoted in Ref. 10). 

Since Scheme III seemed the most applicable method to the transonic regime it was first 
decided to use the collocation method of Jurak et al (Ref, 6) which is similar to Scheme III. By this 
method results were obtained for M as low as 1.1 but even here convergence was slow and the results 
for M < 1,5 did not compare too well with experimental data. Even the Scheme III results given by 
Belotserkovskiy and shown on Figure 6 do not look sufficiently accurate. 

In view of these shortcomings of other techniques it was decided to attempt a relaxation 
solution of the full potential equation for inviscid steady transonic flow. In this equation the flow is 
assumed to be irrotational and hence isentropic. However, the shock is treated as a discontinuity and 
either Rankine Hugoniot jump relations or isentropic jump conditions (neglecting momentum conser- 
vation) are used. Use of the Rankine Hugoniot conditions, applied to flow about spheres, gives shock 
locations which compare reasonably well with experiment at Mach numbers 1.17, 1.30 and 1.62. Ex^ 
periments at lower Mach numbers will be made at NAE in the near future. 

As this paper was in preparation the work of Frank and Zeirep (Ref. 11) came to the authors’ 
attention. In this work they give a semi empirical formula for stand-off distance which is obtained by 
modifying a formula derived for slender bodies of revolution. Generally their prediction of stand-off 
distance is significantly greater than our prediction. 

The first section to follow briefly describes the preliminary attempts made to solve the prob- 
lem. Thereafter the paper is concerned with the potential approach. We first develop, following 
eTameson (Ref. 1), the finite difference scheme and analyse its stability. Then the shock representation 
and shock jumps are considered. Finally the two stage iteration procedure, one for the shock and one 
for relaxation of the interior equations is considered, 

2.0 PRELIMINARY ATTEMPTS 

Several numerical methods for the calculation of the bow shock were tried. They included 
the time dependent approach (Ref. 3), a collpcation method similar to Scheme III of the method of 
integral relations (Ref. 4) and the method of lines as described in Reference 5. 
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In the time dependent method a flow field and shock wave are assumed and these are adjusted 
whilst marching forward in time until a steady state is reached. Good results can be obtained in this 
way for higher Mach numbers but at M less than about 1.3 convergence is slow and the shock wave is 
not smooth. 

The collocation method, used in Reference 6 at M = 10, (this is similar to Scheme III of the 
method of integral relations, Ref. 4) was next attempted. The stand-off distance was found to be rea- 
sonably good for M > 1.5. 

Finally a method of lines solution was attempted and this gave reasonable results for M > 1.3 
although even at this Mach number convergence was slow (Ref. 7). 

Some results of these computations are shown in Figure 6, where we show a logarithm of 
stand-off distance against a logarithm of Mach number. Since these preliminary attempts gave such 
poor results at the lower Mach numbers it was decided to assume the flow irrotational and use the 
potential equation. 


3.0 EQUATION OF MOTION AND THE BOUNDARY CONDITIONS 

3.1 The Region of Computation 

Referring to Figure 1, O is the origin of the spherical polar co-ordinate system (r, 0 ) set back 
a distance d from the leading edge of the body. ODC is the line of symmetry 6 = 0 and OAB is the 
“ cutoff” line 6 = tt/ 2. AD is the axisymmetric body and BC is the stand-off shock wave. Then the 
region of computation is ABCD, 

In Figure 1 we draw a general ellipsoid although we confine ovir attention to spheres in the 
results section. In this case we tsJce unit radius and set d = 1.5, A = 0.5, 

3.2 Transformations 

Let r = G(6 ) and r = F(0 ) be the equations of body and shock respectively. 

In order to fix the shock and body as co-ordinate lines let us introduce new variables (^, a) 

such that 


r - G(6) 
F (6) - G (6) 


0 = aa^ + ba 


( 1 ) 


The region of computation is therefore 0 < ^ < 1, 0 < a < , where corresponds to 0 = 7 t/2 

(the determination of a, b and is shown in Appendix A), 

The above system of co-ordinates could be used as it stands. However one should notice 
that at Mach numbers approaching unity the shock wave is many body radii from the axis (especially 
AB > AO in Fig. 1) and while flow quantities do not change greatly at distances greater than two or 
three radii from the axis they do change considerably on approaching the body. Thus a transformation 
which puts more co-ordinate lines near the body is desirable. Such a transformation is 


^ = A'rj + B'tj3 + C'rj^ r = a 


(2) 


where A', B', C' are chosen as shown in Appendix B. 
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The transformation of the equation of motion which is to follow will be effected by the 

formulas 


where 


9 

9 

9 

9 

9 


” 9| ’ 

90 

Oa + 

® 9a ® 9| 

92 


92 

9 


9r2 

= fc2 

^ 9^2 

’ 3r90 

~ 9| 


92 

9 


9 -> 

— + f " 
9a 

92 

~902 


+ Oee 

9^2 


92 

9|2 


92 

9^9a 


ii 

9o2 


92 

9^9a 


tr = 


F - G 


(G' + ^(F' - G')) 


^re = (F' - rj') 

^ee = -^re (G' + ^(F' - G')) - (G" + ^(F'^ - G") + (F' - G')) 


Oq = (3aa2 + b)'^ 

^ee ~ ~^aal^l 

and by 

A - A _ A 2 ii 

9 ^ ’ 9 t ? 9^2 ■■ ^ 0 r ,2 

where 

77^ = 1/^^ and 


(3) 


(4) 


(5) 


3.3 The Equation of Motion 


In spherical co-ordinates (r, 0) shown in Figure 1, the equation of motion is a2 div 
V = V • grad ('/ 2 V 2 ) where V = (u, v) and a2 is given by 


u 


2 + v2 + 


7-1 


a2 = vi 


7-1 


( 6 ) 
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We assume the flow is isentropic, thus irrotational so that curl V = 0 and hence we can 
define a disturbance potential 0 by V = + V0 where V„ is the free stream velocity vector. 

Making the appropriate substitutions into (6) gives an equation tor <(>: 


( 


1 - 



2uv 

a^r 



^Od 






/ 2uv 


+ cotO 


/ r 


^0 

~2 


as long as 0 0. If 0 


0, (7) becomes 



<^rr 


200 0 

4* 

r2 


20, 

r 


(7) 


( 8 ) 


In these equations is given by 


a2 




(9) 


where 


u = - V„ cos 0 + <pj 


V = sin0 + — 00 
r 


(10) 


Rather than make the substitution of transformation (3, 4, 5) directly into (7) we first con- 
sider the equation written in a form suitable for application of Jameson’s rotated difference scheme 
(Ref. 1). 

For this scheme we need to write the principal part of Equation (7) as 

0SS 0NN 

where S is the local streamline direction at the mesh point and N is normal to the streamline; the sub- 
scripts refer to partial differentiation in the appropriate direction. Now we have by definition 

90 


_9 

9r 


9S 


u — + — 
9r r 


U — + V 
drj 



where 
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U 


= ^ and V = ^ 


(11) 


(N.B. U and V are important physically as their sign indicates the direction of flow in the (rj, r) plane 
as shown in Appendix C), and also 


Therefore 


9 9 u 9 

q T— = -V — + — 


9N 


9r r 90 


- 9 - 9 

= -V — + U — 
9t 7 9r 


where 


V = (vf. - } (,) 


u 


7]t and U = — Og 
^ r 


( 12 ) 


q^ 


92 

= U2 

92 

92 

+ V2 

92 


9t?2 ■" 

drjdT 

9 r 2 

92 


92 — 

92 


92 

9 N 2 

= V2 

9772 " 

drjdr 

+ U2 

9t2 


(13) 


Transforming (7, 8) to the (^, o) co-ordinates gives 


0 ^ 0 : A<p^^ + + D<p^ + E<p„ - 0 

(14) 

/ u2 \ 2 / \ 

^ ^ ^ <^aa + 20g ^ “ j “0 

and transformation to (tj, r) gives 

(l - + I 


+ E0^ = 0 


^r' 


^eo ^ 
r2 r 


(15a) 


(15b) 
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where 


/ 

u2\ ~ 

fc2 

2uv / v2 \ 




- ?) 

" ^ " ^ 2 ) 


(16a) 


2 

/ v2 \ 2uv 



B 

J.2 ^9^0 1 



(16b) 



/ v2 \ y 




C = aj 

(^-rOA^ 


(16c) 


/2uv 

\ ^0 / v2 \ ^60 

2uv ^ 


r 

+ ( + 

V a2 


^27 

a"^r 

(16d) 


E 


( v2 \ / 2uv \ 

1 - - j Oeelr^ + + cot0 j a, /r2 


Now (15a) can be written 


(l - ^) <^>ss + <^NN + +Brt^ (j>^ + E0^ = 0 


(16e) 


(17) 


as shown in Appendix D. Equation (15b), applied only at r - 0, is left in the given form since the 
flow is obviously in the tj direction. 

3.4 Boundary Conditions 

(i) On the body the normal velocity is zero. Thus 


hence 


and so 



( 18 ) 


giving (j>^ on the body (tj = 0). 
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(ii) Line of Symmetry 

In this case we introduce an image line (r = -At) and make the potential symmetric. 

(iii) Along AB 

There is no boundary condition given at the cutoff line so we can only assume that <j> is well 
behaved in this region (no shocks). Then we can introduce a line t = t^ extrapolate 

quadratically to this line from the values of 0 at t^ » ''"m ax “ “ 2 At. 

(iv) Equations at the Shock Wave 

At the shock (17 = ^ = 1) we apply the Rankine Hugoniot (R-H) conditions expressing 
continuity of tangential velocity, energy, mass, and momentum. These are written as 



F' 

F' 

(19a) 

V + 

— u 
F 

= v„ + — u„ 

“ p “ 


pv„ 

= P v„ 

r'oo ’ n 

00 

(19b) 

p + 

pvj 

= P<» + 

(19c)* 

7 P ^ 
-Ip 


Y Poo 

+ V 2 WI 

7'- 1 p„ 

(19d) 


where 


/ 

VI 

These equations are solved for u, v, p and p. In particular u is needed to provide boundary conditions 
for <j) as shown below. 

Firstly we can show that (p is constant along the shock. We know that V*^= -yDecause 
of continuity of tangential velocity. ButV = V„ + V<^andso Oor90/9t = 0. Hence 0 is 

constant along the shock. If we define the constant (which is arbitrary) to be zero then we have 


0 = 0 at T? = 1 


( 20 ) 


or 


We also know the R-H value of u at the shock (u^ say) by solving (19), and hence we have 

Us = +0^ 

= 

Us " U« 

( 21 ) 




We could omit this equation and use p/p'*' = p„ /p„'*' for the isentropic shock relation. 
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4.0 NUMERICAL PROCEDURE 

4.1 Relaxation 

Consider the computational grid shown in Figure 2. Referring to Equations (17) and (15b) 
we write the first derivatives as 






2Arj 


<t>r = 


2Ar 


while second derivatives are written in central difference form if the flow is subsonic (q^ < a^ ) or in 
backward difference form for <^55 if q2 > a^, as follows. 

(i) Subsonic, q2 < a^ 

Following Jameson (Ref, 1) we use 


I'nri 


^ti-1 

Ar?2 


(22a) 


0 


■nr 






i-1 , j 


j+1 *^1+1, j-1 + 


'i+1, j+1 


4Ar?Ar 


(22b) 



for all the second derivatives terms, where co is the over-relaxation factor and is taken as 1.8 in the 
present calculations. The superscript + implies that new values are to be computed (for i subscript) or 
used (for i-1 subscript) in the iteration process, otherwise values of 0 from the previous step are used. 


Thus in the subsonic region we have a tridiagonal system of equations for (j>l^ j , 0^^ j and 
i-i to solved at each step of the line relaxation procedure. 

(ii) Supersonic q2 > a^ 

Considering (1 - q^/a^ ) 0gg + (f>^^ we use (22a) and (22b) for the and 0^^ contribu- 
tions to 0Nf^ while for in the expression we use 


0rr 


<I>1 


i-1 . j 




At^ 
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For contributions to 0gg we use all old values of 0 as follows 


<(>■ 


vn 


AT72 


0 = + 




^17T 


ArjAr 


Ar2 


where the upper sign is used if U > 0 and the lower sign if U < 0. These formulas give a finite dif- 
ference representation consistent with the flow direction since U is proportional to the velocity in the 
r = const direction and V is porportional to the velocity in the rj = const direction (see Appendix C). 

Now to investigate convergence in the same manner as Jameson (Ref, 1) we have to introduce 
an artificial time into the difference equations so that for instance 


i>tj = ^ij + At 

Then in we have 


u (OLD VALUES) y 

— (V0,, - U0,J = 0NN - - 0Nt (23) 


are used, there is no 0gj contribution. Therefore we must artificially 
add some 0gj to the finite difference equation so that we satisfy Jameson’s necessary condition for 
convergence (see Ref, 1 for full details) i,e. 


(OLD VALUES) 

0NN + 

At 


while in 0gg, since all old values 


> |32 (mJ - 1) 


(24) 


where Ml is the local Mach number and -a and jS are the coefficients of 0gj and 0jsj^ respectively; in 
oxir case j3 = -U/q*At/Ar (see (23)), 


U V 

0St " q ^ 


Now 
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which can be written in finite difference form as 


+ H - <^ij + 0i,jTi \ 

q V Ai?At / 

+ X ~ ~ ^ \ 

q \ ArAt ) 


with the upper signs used if U > 0, otherwise lower signs (since U indicates the flow direction). Thus 
the appropriate amount of 0g^ to be added is 


Q 


|U| M 
q At 


^St 


where 


Q2 > - 1 (25) 

to meet condition (24). Since we are dealing with transonic flows, by experience we can expect 

to be about three so that choosing IQI = 4 should be sufficient. As to the sign of Q it should be chosen 

so that diagonal dominance of the tri-diagonal equations is enhanced as follows. 

The coefficient of j_i is 


_ Qu iui 

q2Ai?2 q^ArAT? 

(U > 0) 

the coefficient of <j)*^ j is 


y 2 Q |ui u 

q^Ar? q^ArArj 

(U < 0) 

and the coefficient of ; is 


-2y2 _ Q iui |U| Q |u| y 

q^Arj^ q^Ar^ q^ArA?] q^Ar^ 


Clearly Q < 0 is required since then the magnitude of the diagonal term will be enlarged (unless 
y (= Og v/r) is negative in which case the flow direction is clearly wrong). In practice a test is always 
made to ensure that the system is diagonally dominant; if it is not an error message is printed and com- 
putation ceases. This dominance is necessary to ensure a well behaved solution to the system of alge- 
braic equations. 
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Unfortunately this process does not give convergent results. As an example of the divergence 
we took the converged solution for M = 1.3, P-= 12, N = 7 (i.e. 5 X 10 mesh) and inspected SA0?j 
on successive iterations. On the first iteration SA0?j was 10.59, on the second 67.5 and then no diag- 
onal dominance was present on the next iteration. 

In order to gain some insight into the amount of to be added to give convergence we next 
looked at the amplification factor, G, given by 


G = 




In Appendix E we show that |G,' is less then unity in a simplified case where U = U, V = V and 
\j2 + = q2 , provided 


■ a > q Mf 

Ar ^ 

In practice, however, using a = q At/Ar did not give convergence and we finally had to use 


a = 




+ 


3 



to give convergence. Cases which were very slow to converge were those in which M was close to unity. 
This is probably due to the ‘cutoff’ line, AB in Figure 1, lying partially in the subsonic part of the flow 
field. Table 1 illustrates the convergence when M = 1.02 by listing SA0? j on successive iterations. 
Compare this to the M = 1.3 case shown in Table 2. It can be seen that at M = 1.02 the convergence 
is much slower. Since convergence is slow at the lower Mach numbers fine grids are not computed in 
these cases. 

4.2 Finite Difference Application of the Boundary Conditions 

Referrmg to the boundary conditions in Section 3.4 we now write them in finite difference 

form. 

(i) To apply (18) at 77 = 0 we introduce a line of points at tj = -Atj and write 0^ as 


while 0^ is calculated as 


3 ~ 1 

2Atj 

^i + 1 , 2 “ ^i-1 , 2 

2At 


Substitution into (18) then gives values for 0j j . 
(ii) From symmetry 
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(iii) By quadratic extrapolation 

<^N + l, j “ j - j 

(iv) p = 0. We also need from the shock boundary conditions. Nowwe know 

from (21) thus we could write 




P ‘^’i, P-1 

At? 


and accordingly find 0; p_ j . However, the above formula is only first order accurate for the first 
derivative whilst we really require second order accuracy for consistency with the finite difference 
equations used in the elliptic region. We obtain second order accuracy by using the equation of motion 
(15) to find 0^^ at the shock wave ^0® ^ say^ 

Then we have 


p-i = <^i, p - At?0^ + — 0^^ (26) 

giving, in effect, second order accuracy for the finite difference representation of 0^ at the shock. 

Having seen how to apply the boundary conditions (once a shock shape has been estimated) 
we can now iterate, by relaxation, the equation of motion vioritten in finite difference form. 

4.3 The Inner Loop Iteration 

The preceding describes the finite difference representation of the equations of motion and 
the boundary conditions. To drive the residuals in these equations smaller we solve for j starting 
with line i = 2 (the line of symmetry) and proceeding to line i = N. On each line we solve for (p^ j 
(j = 2 ... P-2 inclusive) by solving the appropriate tridiagonal system of equations. Values for 
<Pi- 1 , j are known except at the line i = 2 and here we are forced to use old values. 

Thus we have a technique for driving residuals smaller at the points i = 2 . . . N, j = 2 ... P-2 
once 0{^ p _2 are fixed. The next section describes how an outer iteration is used to improve the shock 
shape and hence 0j p_j (from (26)) so that the residuals at j = P-1 are made smaller. 

4.4 Newton Iteration of the Shock Wave — The Outer Loop 

We have shown how 0j^ p_j is obtained (26) from the shock wave relations and also we have 
shown how the values of 0 in the rest of the flow field, assuming 0j p_ j is computed from (26), are 
obtained. However we still have to iterate so that the shock wave r = F(0) is adjusted to the hopefully 
correct position. To make the adjustment we must have some outstanding equations still to be satis- 
fied. Clearly these equations are the finite difference form of the equation of motion written at each 
point at the line next to the shock wave, i.e. j = P-1, i = 2 . . . N. Unfortunately, these equations 
cannot be written in backward difference form in the supersonic region since we only have values 
available at one mesh line further upstream — thus a central difference form has to be used even in the 
supersonic region. This, however, will not affect convergence since a Newton method is to be used on 
this line to drive the equation residuals, Rj (i = 2 . . . N), smaller by adjusting the parameters which 
define the shock wave. These parameters Aj (j = 0, 1 . . . NSH-1) are defined in the next subsection. 

In the outer iteration the parameters Aj are changed by an amount X6Aj where 5Aj is given by 
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NSH-1 

( N 9% 

I 

N 

2 

j = 0 

2 

[ k = 2 3Aj 

3A| ) 

i 


Rfc 


3% 

3Ai 


i = 0, 1 . . .NSH-1 


'y 

The above is the Newton step to minimize Sj = 2 Rj with respect to A; . X (0 < X < 1 ) is 

i = 2 

chosen so that Sj decreases (normally X = 1 but may have to be smaller to ensure a decrease in Sj ). 

The partial derivatives in the above equation are found by differences using a perturbation to Aj of 10“'* * 

4.4.1 Form of the Shock Wave 


We could let the shock wave r = F(0) be defined by its values F (02 )> F(^ 3 ) . . . F(0f^ ). 
However there are two reasons against this form. This first reason is for efficiency — if we can define 
the shock by less than N parameters then the Newton iteration of the shock is computationally quicker. 
More importantly, though, the finite differences approximations to determine F' (0) and F” (d) will be 
poor in the region near 6 = 90° because, here, there are large gradients of F since F is going rapidly to 
infinity at the lower Mach numbers. Thus it is important to have an analytic form for the shock that 
can be differentiated analytically. 

Such a form is given by 


r = F(0) = 


1 - cosd' 
cosO - cosd' 


NSH-1 

2 

0 


A„« 


2n 


where 6' is the complement of the Mach angle. This has the correct form since such a shock wave is 
asymptotic to the Mach line and is also symmetrical as required by the problem. Note that F (0) gives 
the stand-off distance except for a constant which depends on the location of the origin of the co- 
ordinate system. In the cases computed, this origin is chosen such that A = VzDE (see Fig. 1) where 
E is the centre of the sphere. 

4.5 Iteration Procedure 

In the iteration procedure described previously the inner loop drives S 2 = 2A(^j j smaller 
by relaxation while the outer loop adjusts the shock shape so that residuals Rj at the line next to the 
shock are driven smaller. One problem in this scheme is knowing when to adjust the shock shape — 
if we adjust it when 83 is still large then the adjustment may be completely meaningless while on the 
other hand if we wait until 83 is too small the efficiency of the scheme suffers. Thus it was decided 
that the shock would not be adjusted until 82 was less than P*N*10“'* . Also to avoid possibly meaning- 
less changes to shock position a change in F (0) of greater than 0.05 on each step was not allowed. 
Likewise for efficiency a change in F (0) of less than 0.001 was not made. A typical iteration is shown 
in Table 3 at M = 1.3, mesh size 5 X 10, i.e. N = 7 and P = 12, and 0jj set to zero initially. 

A flow diagram of the iteration scheme is shown in Figure 3. 


5.0 80ME PROGRAM DETAILS 

5.1 Initial Estimates 

We first selected M = 1.3 as a starting Mach number, made a rough estimate of the coeffi- 
cients that form the shock, set j = 0 on a 5 X 10 mesh and let the process run. Having obtained 
this solution we then used the results to give the estimates for the Mach numbers 1.17 and 1.10. 
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Having obtained these solutions we then extrapolated quadratically the shock coefficients and the 
initial down to lower Mach numbers. In this way solutions were obtained for M = 1.08, 1.06, 1.05, 
1.04, 1.03 and 1.02. At some of these Mach numbers finer grids (10 X 20 and 20 X 40) were computed. 
To get estimates for the finer grids we used the same shock values from the coarser grid and interpolated 
0jj linearly to obtain values at the extra mesh points. 

5.2 Computer Time 

The three meshes 5 X 10, 10 X 20 and 20 X 40 took approximately six minutes on an IBM 
360/67 for the M = 1.3 solution. At M = 1.05, five minutes were required to compute the coarse 
and medium grids. 


6.0 RESULTS FOR FLOW ABOUT A SPHERE 

6.1 Program Checks 

As a check on the program a mass balance was made to see how the mass flow across the 
circle of radius HK (Fig. 1) compared with that across AB. Table 4 shows the total mass inflow to- 
gether with the mass outflow. Good agreement is noted. 

As a further check comparisons were made with the results of South’s axisymmetric program 
RAXBOD (Ref. 2) which treats the whole flow field extending to infinity and does not allow for dis- 
crete shocks jumps; ‘shocks’ appear as discontinuities in the velocity components. Results of the com- 
parisons are shown in Figures 4 and 5. These show Mach number distributions along the body and also 
along the axis of symmetry for free stream Mach numbers 1.3 and 1.02. The M = 1.02 results have not 
converged to the preset accuracy requirement. However it can be seen that agreement is good when 
fine meshes are used with the exception of Mach number behaviour in the supersonic zone along the 
body. Here first order accuracy is used so that a larger discrepancy is expected and also the cutoff Line 
(not used in RAXBOD) may have an effect on the SHOCKFIT results. 

6.2 Residuals 

At the end of each computation one extra pass is made to calculate residuals at each point 
(with factor a set to zero). The size of the residual does not mean too much unless we compare it 
with a relevant quantity. Here we compare residuals with the right hand side (d) of the equation 


a(/)^_l + + c0^ + i 


= d 


which is the finite difference for^i of the equation of motion. The residuals and right hand sides are 
printed in Tables 5 and 6 for K^ach numbers 1.3 and 1.02 (the latter at the fine mesh have not con- 
verged). Note that the residuals at the line next to the shock are much smaller at M = 1.02 than at 
M = 1.3. This is probably due to the entropy change across the shock being much bigger at the larger 
Mach number and this change is inconsistent with the isentropic assumption for the rest of the flow 
field. 

6.3 Shock Wave Location 

The main results of the paper are given in Table 7 which shows the distance F (0) of the 
shock wave from the co-ordinate origin for polar angles 6 of 0, 18, 36, 54, 72 and 90°. It can be seen 
that results depend quite strongly on mesh size. The main problem is to be sure of accuracy at 0 = 0° 
since it is here that the correlation with experimental results will be made. 

First let us see how much accuracy is needed in F (0) so that the free stream Mach number 
can be determined to an accuracy AM. To do this we need to know an approximate relation between 
F (0) and M. Such a relation is given in NACA 2000 (Ref. 8) i.e. 
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F (0) = - (M- + 0.5 

3 

(see Fig. 6 to compare this form with the numerical computations). A change of less than AM in M 
thus requires 


AF (0) < - (M- AM . 

Taking AM = 0.002 we can deduce limiting values for AF (0) as shown in Table 8. The last column 
in the table shows the limiting value for Ain (F (0) - 0.5) and can be calculated as 

1 AM 

Ain (F (0)- 0.5) . 

^ 3 M- 1 

The size of this latter quantity is marked by a symbol I on Figure 6 which shows In (F (0) - 0.5) 
versus - In (M - 1). 

Now observing Table 7 we see that the 5 X 10 and 10 X 20 results for stand-off distance 
(3 F (0) - 1.5) differ by quite a large amount so that an extrapolated zero mesh size solution should 
be used. This latter solution is obtained by extrapolating the 5X10 and 10 X 20 solutions by assum- 
ing the error is of order mesh size squared — the resulting values are given in Table 7 by the name 
EXTl. It can be noted that at M = 1.3 the EXTl result is practically identical to that obtained by 
carrying out a similar extrapolation from the 10 X 20 and 20 X 40 grids (called EXT2 in Table 7). 

At M = 1.06 the difference in F (0) between the EXTl and EXT2 results is 0.01 — well within the 
accuracy requirement given in Table 8. Thus we can have some confidence in the accuracy of extra- 
polation from the 5X10 and 10 X 20 results and so we propose that the EXTl results are our final 
results with accuracy sufficient to be used to correlate stand-off distance to Mach number. We use 
this procedure rather than compute further solutions on the 20 X 40 grid since computation is rather 
expensive on the fine grid ($88 at M = 1.06). 

Using the EXTl results we note that, on the log scales of Figure 6, they lie almost on a 
straight line given by 

In (F (0)- 0.5) = -0.505 In (M- 1) - 0.045 

or F(0) = 0.5 + 0.956 (M- 1)"°'^°^ (27) 


The values given by this formula are also tabulated in Table 7. Notice that the differences between 
F (0) given by (27) and F (0) given by EXTl, for M < 1.08, are much less than the acceptable toler- 
ances given in Table 8 — given some confidence in our straight line fit. 

The RAXBOD results for stand-off distance, given in Figure 6, are found by interpolating 
between Mach numbers calculated along the axis of symmetry to find at what distance from the body 
the Rankine Hugoniot value of Mach number is realized. This procedure should give reasonable accu- 
racy since Mach numbers calculated at some distance from the calculated shock discontinuity are con- 
sidered accurate and the discontinuity, from experience, is always upstream of the real shock jump. 
Plots of M(axis) from RAXBOD are shown in Figures 4b and 5b. 

In the work of Frank and Zeirep (Ref. 11), who use a modification of a formula derived for 
slender bodies, is given the formula 

/ 0.14 (7 + 1) m2 
F(0)- 0.5^2 ( 
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Values obtained by this formula are plotted on Figvure 6. It can be seen that the slope of this curve is 
significantly greater than the slope of our curve thus predicting much larger stand-off distances when 
M < 1.10. Hopefully experiments at these low Mach numbers, to be performed shortly at NAE, will 
settle the discrepancy. 

The shock shapes are compared with experimental data from NACA 2000 (Ref. 8) in Figure 7. 
It can be seen that satisfactory agreement is obtained. 


7.0 CONCLUSIONS 


We have used flameson’s rotated difference scheme combined with a Newton iteration of the 
shock wave to obtain a prediction of the shock wave stand-off distance from a sphere. The same method 
could also be applied to more general axisymmetric bodies. 

As far as is known this is the only rigorous theoretical work which predicts the flow field 
solution and shock wave location for Mach numbers less than 1.05. Frank and Zeirep in their predic- 
tion use a semi-empirical formula based on a modification of slender body theory while Hsieh, using 
the time dependent approach, takes more than 22 hours for a (nonconverged) solution at M = 1.05. 

The present method needs about five minutes at this Mach number. 


The accuracy of the present method can be seen to be fairly good at Mach numbers 1.62, 
1.3 and 1.17 at which results have been compared with experiment. Experiments to be performed 
shortly at NAE at the lower Mach numbers will be used for a further comparison. 

Our results (referring to stand-off distance) of best accuracy, given by extrapolating to zero 
mesh size from a 5 X 10 and a 10 X 20 grid, lie practically on a straight line when plotted on the 
log scales of Figure 6. The resulting formula for F (0) is given by (27). 
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TABLE 1 


VALUES OF S 2 AND Aq DURING THE ITERATION PROCESS M = 1.02 


ITN# 

5 X 10 Mesh 

10 X 20 Mesh 

20 X 40 Mesh 

S 2 

Ao* 

S 2 

Ao 

S 2 

Aq 

0 

1.6, -7 

6.553 

3.6, -3 

6.666 

1.3, -3 

6.941 

5 

3.6, -6 


2.0, -2 


1.8, -3 


10 

4.6, -7 


1.1, -3 


2.1, -4 


15 

1.2, -6 


3.4, -4 


7.3, -5 

6.995 

20 

8.8, -7 


1.0, -4 


7.7, -5 


25 

4.9, -7 


4.5, -5 


8.6, -6 

7.017 

30 

2.5, -7 


2.4, -5 


2.9. -5 


35 

1.3, -5 

6.532 

1.3, -5 

6.722 

1.0, -5 


40 

3.3, -6 


3.4, -5 


1.0, -5 

7.029 

45 

1.1, -6 


3.7, -6 

6.747 

1.1, -5 


50 

5.2, -7 


9.1, -6 


Iteration terminated 

55 

2.6, -7 


1.1, -5 

6.791 



60 

1.3, -7 

6.582 

1.5, -5 




65 

3.3, -6 


1.6, -5 

6.841 



70 

2.8, -7 


1.8, -5 




75 

2.0, -6 

6.632 

1.6, -5 

6.891 



80 

9.2, -7 

•> 

1.7,~5 




85 

1.6, -6 

6.661 

1.6, -5 

6.941 



90 

2.4, -7 


1.7,-5 




95 

6.2, -8 

6.666 

Iteration terminated 



100 

3.3, -8 







* Shock wav'€ coefficients were changed during the five iterations previous to the listed value. 








PJ^KcmrNCr 
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TABLE 2 

VALUES OF SA0?j AND Aq DURING THE ITERATION PROCESS M = 1.3 


ITN# 

5 X 10 Mesh 

10 X 20 Mesh 

20 X 40 Mesh 


A„ * 
•^0 

S 2 

■^0 


Aq 

0 

2.3, -1 

2.5 

1.9,-3 

2.397 

4.4, -4 

2.414 

5 

1.7, -1 


9.4, -5 


7.6, -5 

2.416 

10 

3.4, -2 


9.3, -5 

2.408 

1.1, -5 

2.419 

15 

5.5, -4 


8 . 6, -6 

2.411 

6 . 2, -6 


20 

2.6, -5 


5.8 , -6 

2.413 

2 . 1, -6 


25 

6.3, -4 

2.450 

2 . 2, -6 


1.3,-6 


30 

3.5, -6 

2.411 

2 . 0, -6 

2.414 

6.3, -7 


35 

8.3 , -6 


1.8, -7 


3.6, -7 


40 

1.5, -5 

2.402 





45 

2 . 1, -6 

2.398 





50 

2.2, -7 

2.397 





55 

2.5 , -8 







* Shock wave coefficients were changed during the five iterations previous to the listed value. 


TABLE 3 


PROGRESS OF ITERATIONS M = 1.3, 5X10 MESH, (^i j = 0 INITIALLY 


Shock Change 
Made on ITN # 



Shock Coefficients 

S2 

Si 

^0 

Ai 

^2 

■1 

A 4 

As 

0 

2.3, -1 


2.5 

- 0.1 

0 

0 

0 

0 

24 

4.6 , -6 

1 . 0, -2 

2.450 

-0.1184 

0.0514 

-0,0270 

0.0082 

- 0.0010 

30 

3.5 , -6 

1.3, - 3 

2.411 

-0.1251 

0.0773 

-0.0370 

0.0122 

-0.0018 

37 

2 . 8, -6 

2.0, -4 

2.402 

-0.1235 

0.0785 

-0.0377 

0.0138 

- 0.0021 

42 

1 . 8,-6 

9.9, -5 

2.398 

-0.1242 

0.0810 

-0.0394 

0.0145 

- 0.0021 

46 

9.2, -7 

7.9, -5 

2.397 

-0.1242 

0.0811 

-0.0391 

0.0142 

- 0.0020 

55 

2.5 , -8 


No Further Change 
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TABLE 4 


CHECK ON MASS CONSERVATION 


M 

Mesh 

Mass Flow 

Across HK 

Across AB 

1,62 

5 X 10 

7.90 

8.04 

1.3 

5 X 10 

23.10 

23.09 


10 X 20 

23.24 

23.26 


20 X 40 

23.28 

23.31 

1.7 

5 X 10 

56.74 

56.71 


10 X 20 

57.93 

57.93 

1.10 

5 X 10 

136.3 

136.2 


10 X 20 

144.7 

144.7 

1.06 

5 X 10 

335.1 

334.9 


10 X 20 

345.3 

345.2 

1.04 

5 X 10 

710.9 

710.7 


10 X 20 

597.1 

696.9 

1.02 

5 X 10 

2774 

2773 


10 X 20 

2659 

2659 


20 X 40 

2660 

2659 (not converged) 











TABLE 5a 

RESIDUALS AND RIGHT HAND SIDE OF FINITE DIFFERENCE EQUATIONS 
AFTER ITERATIONS COMPLETED M = 1.3, 5X10 MESH 



2.5,-6 

7.8,-l 



8 

-2.2, 

-6 

1.6 


-5.2, 

-7 

7.3, 

-'1 

-1.7, 

-6 

5.8, 

-1 

-5.6, 

-8 

2.4, 

-1 

9.1, 

-5 

4.1, 

-2 

-5.0, 

-4 

1.3, 

-2 



(Line Next 
to Shock) 

10 

11 

-2.9, -6 
2.0 

-2.5, -3 
-2.7, -1 

-1.3, -6 
1.5 

1.7,-3 
-3.1, -1 

-1.8,-6 
9.7, -1 

8.9, -4 
-1.8, -1 

-7.3, -6 
5.0, -1 

3.1, -5 
-9.7,-2 

9.3, -5 

2.3, -l 

-2.8, -3 
-5.1, -2 

-5.7,-4 

1.2,-1 

-9.4, -3 
-1.5,-2 


















TABLE 5b 


RESIDUALS AND RIGHT HAND SIDES OF FINITE DIFFERENCE EQUATIONS 
AFTER ITERATIONS COMPLETED M = 1.3, 20 X 40 MESH 


Line# 


(Body) 
J = 2 

6 

10 

14 

18 

22 

26 

30 

34 

38 

(Line Next 
to Shock) 

41 

2 

RES 

RHS 

-3.9, -4 
4.0, 2 

-3.7,-4 
1.0, 2 

-3.5,-4 
9.1, 1 

-3.2, -4 
7.4, 1 

-2.9, -4 
5.7, 1 

-2.5, -4 
4.1, 1 

-2.1, -4 
2.8, 1 

-1.8, -4 
1.7, 1 

-1.5,-4 

8.7 

-l.l,-4 

3.1 

-1.5, -2 
8.6, -2 

6 


-4.1, -4 
2.8, 2 

-6.1, -5 
4.6, 1 

-5.8,-5 
4.2, 1 

-5.3,-5 
3.4, 1 

-4.7,-5 
2.6, 1 

-4.1, -5 
1.9, 1 

-3.5, -5 
1.3, 1 

-3.2, -5 
7.6 

-3.6, -5 
3.9 

-4.9, -5 
1.2 

-1.2, -2 
-2.1,-1 

10 


-1.1, -3 
1.3, 2 

-9.1, -5 
1.7, 1 

-8.9,-5 
2.1, 1 

-8.4, -5 
2.0, 1 

-7.8, -5 
1.8, 1 

-7.3,-5 
1.4, 1 

-6.9, -5 
9.3 

-6.7,-5 

5.6 

-6.9, -5 
2.8 

-7.3,-5 
8.7, -1 

-1.2,-2 

-1.2,-1 

14 


-2.4, -3 
3.6, 1 

-3.5, -4 
1.0, 1 

-2.4,-4 

7.1 

-1.6, -4 
4.9 

-l.l,-4 

3.5 

-9.3,-5 

4.2 

-8.7, -5 
3.6 

-9.2, -5 
2.5 

-l.l,-4 

1.1 

-1.5,-4 
3.2, -1 

-6.4, -3 
-9.5, -2 

18 


-5.9, -3 
-1.6, 1 

-1.5, -3 
9.2 

-1.1, -3 
3.9 

-8.9, -4 
1.7 

-7.3,-4 

7.7,-l 

-5.6,-4 
4.0, -1 

-3.5, -4 
2.3, -1 

-2.0, -4 
1.2, -1 

-2.1, -4 
3.8, -2 

-2.9, -4 
-3.3,-2 

-1.6, -2 
-8.7, -2 


22 


-1.6, -2 -5.6, -3 -3.8,-3 -2.8,-3 -2.2,-3 -2.0, -3 -1.9, -3 -1.8,-3 

-2.9, 1 3.9 6.6,-l 3.3,-l 2.5, -1 1.5, -1 6.4,-2 1.7,-2 


-1.5, -3 
-3.7, -2 


-3.7,-2 
-3.8, -2 
























TABLE 6a 


RESIDUALS AND RIGHT HAND SIDES OF FINITE DIFFERENCE EQUATIONS 
AFTER ITERATIONS COMPLETED M = 1.02, 5 X 10 MESH 


Line# 


1 

(Body) 

J = 2 

3 

4 

5 

1 

6 

7 

1 

8 

9 

10 

(Line Next 
to Shock) 

11 

2 

RES 

-1.4, -4 

-5.8, -4 

-2.2, -4 

-9.6, -5 

-4.4, -5 

-1.9, -6 

3.2, -5 

3.1,-5 

8.2, -6 

1.2, -6 


RHS 

1.4, 2 

2.6, 1 

8.3 

4.1 

2.0 

9.4, -1 

4.1, -1 

1.7, -1 

6.7, -2 

1.8, -2 

3 


5.8, -5 
1.1, 2 

-2.7, -4 
2.0, 1 

-3.0, -5 
4.8 

1.3,-5 

2.0 

1.6,-6 
9.4, -1 

-6.3, -6 
4.2, -1 

6.1, -6 
1.8, -1 

1.7, -5 
7.1, -2 

5.4, -6 
2.9, -2 

1.1, -6 
6.9, -3 

4 


1.1, -4 
6.3, 1 

-1.5,-4 
1.2, 1 

-1.8, -5 
3.5 

7.6, -6 
1.6 

2.4, -6 
7.3, -1 

-2.1, -7 
3.0, -1 

8.0, -6 
1.2, -1 

1.8, -5 
4.3, -2 

4.3, -6 
1.8, -2 

1.0, -6 
3.5, -3 

5 


1.2, -4 

1.2, 1 

-4.6, -5 
4.3 

7.1, -6 
1.2 

6.9, -6 
8.2, -1 

7.9, -7 
3.6, -1 

4.2, -6 

1.3, -1 

1.3, -5 
4.5, -2 

2.2, -5 
1.5, -2 

1.1, -6 
6.8, -3 

1.6,-6 
6.3, -4 

6 



4.2, -5 
8.5, -1 

-1.9, -5 
4.0, -1 

-1.9, -5 
7.7, -2 

-6.0, -6 
5.9, -2 

9.9, -6 

1.9, -2 

2.6, -5 
5.0, -3 

3.2, -5 

1.2, -3 

-6.6, -6 
1.4, -3 

-1.9, -6 
-3.1, -4 

7 



3.4, -4 
3.7, -1 

5.2, -5 
1.8, -2 

9.4, -6 
2.7, -4 

1.6, -5 
-9.5, -4 

3.4, -5 
-4.9, -4 

3.2, -5 
-2.2, -4 


2.1, -5 
-1.2,-4 

3.8, -6 
1.4, -4 

-7.5, -7 
-7.3,-5 



TABLE 6b 


RESIDUALS AND RIGHT HAND SIDES OF FINITE DIFFERENCE EQUATIONS 
AFTER ITERATIONS COMPLETED M = 1.02, 20 X 40 MESH (NOT CONVERGED) 


Line # 


(Body) 
J = 2 

6 

10 

14 

18 

22 

26 

30 

34 

38 

(Line Next 
to Shock) 

41 

2 

RES 

RHS 

-3.1, -2 
6.6, 2 

-3.0, -2 
1.8, 2 

-2.6,-2 
1.2, 2 

-2.0, -2 
7.6, 1 

-1.1, -2 
3.9, 1 

-3.7, -3 
1.7, 1 

-2.4, -3 
7.1 

-9.6, -4 
2.8 

5.0,-4 

1.0 

1.4, -4 
3.0, -1 

-2.2, -5 
5.1, -2 

6 


-2.3,-2 
4.9, 2 

-6.7, -4 
8.4, 1 

-6.4, -4 
5.7, 1 

-4.9, -4 
3.6, 1 

-5.3,-4 
1.8, 1 

-2.9, -4 
7.8 

-2.7, -4 
3.2 

-5.8,-4 

1.2 

9.3, -5 

4.4, -1 

9.7, -5 
1.3, -1 

2.0, -6 
2.1, -2 

10 


-3.6, -2 
2.4, 2 

-1.9, -4 
3.8, 1 

-4.0, -4 
3.6, 1 

-4.1, -4 
2.6, 1 

-3.7, -4 
1.3, 1 

-1.9,-4 

5.5 

-2.4, -4 
2.1 

-5.2, -4 
7.5,-l 

6.5, -5 

2.6, -1 

6.1, -5 
7.3,-2 

-2.1, -5 
1.1, -2 

14 


-4.0, -2 
3.8, 1 

2.5, -4 

1.5, 1 

1.0, -4 
8.0 

-2.4, -4 
1.1, 1 

-3.2, -4 
6.4 

-1.3, -4 
2.4 

-2.3,-4 
8.1, -1 

-3.9, -4 
2.7, -1 

5.2, -5 
8.7, -2 

-1.1, -5 
2.3, -2 

2.1, -6 
3.0, -3 

18 


-3.3, -2 
-5.7, 1 

-4.8, -5 
1.2, 1 

-7.7, -5 
2.5 

-8.5, -5 
7.2, -1 

-2.0, -5 
1.2 

-2.6, -4 
3.9, -1 


-1.7,-4 
3.2, -2 

-7.6, -6 
9.3, -3 

-3.4, -5 
2.0, -3 

-2.3, -5 
-3.9, -5 

22 


-1.4, -2 
-2.4, 1 

-1.2,-3 

7.9,-l 

-1.1, -3 
8.0, -2 

-2.7, -4 
1.3, -3 

1.7, -4 
2.2, -2 

-2.7, -4 
3.9, -3 


-7.3, -5 
-2.2, -6 

-5.5, -5 
-7.7,-5 

-2.7, -6 
-8.7, -5 

-7.5, -5 
-7.0,-5 
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TABLE 7 

DISTANCE OF SHOCK FROM ORIGIN FOR DIFFERENT d VALUES, 
i.e.F(e). EXT1,EXT2 = VALUES EXTRAPOLATED 
FROM 5 X 10, 10 X 20 AND 10 X 20, 20 X 40 GRIDS 



I 


5 X 10 

5 X 10 
10 X 20 
20 X 40 
EXTl 
EXT2 

5 X 10 
10 X 20 
EXTl 

5 X 10 
10 X 20 
EXTl 

5 X 10 
10 X 20 
EXTl 

5 X 10 
10 X 20 
20 X 40 
EXTl 
EXT2 

5 X 10 
10 X 20 
EXTl 

5 X 10 
10 X 20 
EXTl 

5 X 10 
10 X 20 
EXTl 

5 X 10 
10 X 20 
EXTl 


F(0) by 
Formula (27) 


10.283 
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TABLE 8 

REQUIRED ACCURACY IN SHOCK LOCATION 
TO GIVE MACH NUMBER ACCURATE TO AM = 0.002 


— 

M 

Accuracy Required 
in F(0) 

Accuracy Required 
in In (F(0)- 0.5) 

1.3 

0.004 

0.002 

1.10 

0.019 

0.007 

1.08 

0.026 

0.008 

1.06 

0.038 

0.011 

1.06 

0.048 

0.013 

1.04 

0.065 

0.016 

1.03 

0.095 

0.022 

1.02 

0.163 

0.033 

1.01 

0.412 

0.066 








FIG. I'- CO-ORDINATE SYSTEM 
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i = 2 3 4 N 


FIG. 2: THE COMPUTATIONAL GRID 
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SELECT N,P TO FiX GRID SIZE, E. G. N=7, P = 12 IMPLIES A 5x10 MESH 

SET ITMAX = MAX. NO. OF ITERATIONS 

ESTIMATE SHOCK COEFFICIENTS A^, A^ ••• A^^^^ I 

SET S = P. N. 10'^, SFINAL = P.N. I0‘^ 

ITN = 0 




PROGRAMME 


FIG. 3: FLOW DIAGRAM SHOWING THE INNER AND OUTER 

ITERATION LOOPS 
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4b. MACH NUMBER VERSUS DISTANCE ALONG AXIS OF SYMMETRY 


FIG. 4 : 


RAXBOD COMPARED WITH PRESENT RESULTS, M=l-3 
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FIG. 5: RAXBOD COMPARED WITH PRESENT RESULTS, M= 102 


EXPERIMENT 

f BELOTSERKOVSKIY 


NACA 2000 


.o _ ACCURACY 
REQUIRED TO 
GIVE AM<0 002 



F ( 0 )=%( M -|)'*'»-'/2 
(NACA 2000 FORMULA) 


(l'06) I, 


THEORY 

A 5 X 10 7 
^ 10 X 20 > SHOCKFIT 
0 20 X 40J 

0 EXTRAPOLATED FROM 5 X 10 
AND 10 X 20 MESH 


♦ □ 


CT RAXBOD 36X36 


□ BELOTSERKOVSKIY - SCHEME 


X FRANK AND ZEIREP 


0 ' I ' ' 2 ‘ ' 3 ' '4 ' 5 -In (M-|) 

1-62 1-3 M7 HO 106 104 102 101 — ► M 


FIG. 6: SHOCK WAVE STAND-OFF DISTANCE VERSUS MACH NUMBER 




FIG. 7; COMPARISON OF SHOCK SHAPES 


- 37 - 


APPENDIX A 


The 6-^0 Transformation 


Consider part of the transformation (1) i.e. 


d = ao^ + ba . 


Let 


a = at 0 


7T 

2 


hence 


^ = aa^ + boM 


(Al) 


Now we may require to bunch up 0 = const co-ordinate lines near the axis of symmetry 0=0 (this 
is done in the case of ellipsoid calculations). To do this let 


where a is a constant fixed by the user to bunch lines closer together, thus 


3aaJi 


+ b = ab 


We fix b = 0.1 and solve (Al) and (A2) for a and 0 ]yi giving 




37T/2 
(2 + a) b 


and 


(a - 1) b 
a = 


(A2) 


(A3) 


(A4) 




i'Aun. bLAiNK NUi' I'lLiVlED 
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APPENDIX B 


The ^ T? Transformation (2) 


We require 


t7 = 1 at ^ = 1 

and also, to bunch the qo-ordinate lines near the body, 

ei 

where gj and g2 determine the density of the co-ordinate lines near the body and shock relative to the 
density in the middle. Since the shock departs rapidly from the body as M approaches unity it was 
felt that g] should accordingly grow smaller with Mach number decreasing while g2 is fixed at 0.8 so 
that ^ co-ordinate lines are not too sparse near the shock wave. Thus we choose 



gj = M - 1 

and g2 =0.8 

Now = A' + 3B'7?2 -h 4 C'tj 3 

and so the equations for A', B', C' are 

A' + B' + C' = 1 
A = gj ^ A + 3Bt7^y2 

A' + 3B' + 4C' = g2 ^A' + 3 B'tj {/2 +4 C'i?J/ 2^ " ~ 

where is the value of 17 at | = 1/2 hence 

1/2 = A'r]Yj2 ^^ 1/2 


17 j j2 is found using Newton’s method while A', B', C' are found from linear analysis. 

PRECEDING PAGE BLANK NOT FILMED 
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APPENDIX C 


Velocity Components in the (17, t) Plane 


As defined previously (u, v) are the velocity components in the (r, d) directions. Now let u' 
be the velocity along r = const and v' be the velocity along tj = const. Thus u', v' are the velocity 
components along the mesh lines in the (77, t ) plane and as such indicate which backward finite differ- 
ence expressions to use in the supersonic region. 

To find u', v' we let ju be the angle between the line 77 = const and the line 0 = const. Then 


and 


u = u' + v^ cos /i 
V = v' sin p 


Now the line 77 = const is given by 


and is therefore given by 


0 = 577 = 77 j 5^ = + 8 q 89 ) 

5x _ ^0 

~~ 


so that 


tan/i = 


x8d 


5 r 


Thus 


v' = 


and 


u' = u-v' cosju = u-v COtjU = 




e 


u + V — 



referring to ( 11 ) it can be seen that u', v' are related to U, V and in fact have the same sign since 77^ 

(by definition) and (= 1 /(F - G)) are positive. Thus U and V are used to determine the direction of 
flow in the supersonic region. 
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Consider 


Now 


and 


Hence 


APPENDIX D 


Transformation of the Principal Part to Streamline Co-ordinates (S, N) 


(l - ^)0ss «NN =(A - ^) (u20,„+2UV«„+V2«„) 


i (v20„ - 2UV^„ + 


( 


\J2 + V2 \J2 


q2 a2 / 


\ /uv-uv uv\ 

) ^ ( q2 " a2 ) 


/ y2 + U2 _ V2\ 

\ q2 a2 / ^ 


2 g2 f 


2 / ^nr 


U2 + V2 =(uf, + i,y nl +(v{, - j- ?,) 




= 77^ (u2 + v2 ) |2 + 


u2 + v2 




= TjJq2 


i<!4) 


U2 = (u2|2 + f ^2 ^ 2UV|, 


U2 + V2 U2 


^2 (. 2,^1 ^ y 2 _±_ .2 ^ 

■^r2 a2 a2r2 a2 r / 


- «2 
- 


[>; (■ ■ 3 


2uv 




r r'^ 


(■-31 


At?; 


PBSJCBDING P^GB 


not filmed 


(see (16a)) 
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The 20^^ coefficient is 



Finally 


and so 



(see (16b)) 



(see (16c)) 



^SS *^NN ~ 




+ C<l>rr 


and formula (17) is recovered. 



APPENDIX E 


Evaluating the Amplification Factor 


The equation of motion with the term added is 


( 


1 “ 7 1 ^SS ^NN “ L.O.T. 


(lower order terms) 


(El) 


where 

q^<fss = + V2^„ 

q"*NN - ^ 

For the difference equations in we use 

= (pli + l - 24>ij + <Pii - 1 

4Ar?Ar<^^^ = j , j _ j ~ <Pt-i,i + i " + i , j - 3 + i , j + i 

Ar20^^ = - <Pti - 0ij + <^i + i,j 

while in 0gs we use 

At? 20^^ = 0ij - 20ij _ 1 + 0i j - 2 

At7At0„i = j - j - 1 - - 1 , j + - 1 , j - 1 

Ar20^^ = 0i, j - 20i _ ij + 0i _ 2,j 

assuming, for the analysis, that U > 0 and V > 0. 0§, is represented as 

U V 

<^St ~ ~ ^Vt ~ ^rt 

_ U / - 3 - <^i j + </>i j - 3 \ 


(E2) 

(E3) 


(E4) 

(E5) 

(E6) 

(E7) 

(E8) 

(E9) 


ArjAt 


V / (t>t] - 3 , j - <^i j + 0i - 3 , j 

— 


ArAt 
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In order to analyse the convergence of the iteration scheme the above equations have to be 
simplified. This is done by assuming that U = U, V = V and = q2 . An insight into this 

simpler case may lead to convergence in the more general case. 

Let 0ij = (ElO) 

where G is the amplification factor given by 




A convergent scheme requires |G1 < 1. Now let 

IX = mAi? and v = nAr (Ell) 

and substitute the forms (ElO), (Ell) into (E4) - (E9) and substitute those into (El) to give 


(l - + (l- 


e-‘M _ e' ‘ 


i.) 


2UV 

At?At 


V2 




{( e ‘'^-2 + e ">^) 


GV2 

Atj2 




2UV 

4AtjAt 


- j (g - Ge"**^ - l + e“*^) 


U 


AtjAt 


where 


+ (g - Ge"’"' - 1 + e"'" ) ^ j = q2 • L.O.T. 

/ \ 

\qAt j 


a, = q2 


What is the form to be used for ceAr/qAt? Condition (25) indicates that we use 

cxAt _ Q[U| 
qAt q2 


(E12) 
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where 


1; therefore Q = is taken for the analysis and so 


^ _ Ml lUl ^ 0.1 
qAt q2 


(E13) 


is sufficient. 

Ignoring the L.O.T., since they vanish as At 7 , At 0 using At] /A r = 1, and taking various 
values for Ml > V, /x and v, |G| is then calculated from (E12). The computer program to calculate 
IGl is shown at the end of this appendix. Using the form (E13) it was found that |G| was less than unity 
when Ml = 1.001 but at Ml = 1.834 and higher Mach numbers many of the |G1 values were greater 
than unity. Thus an additional 0§t which is proportional to M^ is indicated and so we next used 



(E14) 


and found |G| < 1 in all cases. 


\ 

On using the form (E14) in the relaxation program it was found that all cases did not con- 
verge — presumably due to boundary conditions or to the simplified model used in the above analysis. 
In fact in the relaxation program it was found necessary to use 


uAr 

qAt 



to ensure convergence in the majority of cases. However, for values of M close to 1 and with a fine 
grid, convergence was sometimes not achieved to sufficient accuracy. For example Table 1 shows the 
convergence history for M = 1.02. The lack of convergence in the latter cases is probably due to the 
‘cutoff’ line AB (Fig. 1) lying partially in the subsonic part of the flow field. 

It may be of interest to note that 

^ |U| + 0.1 

qAt q2 


gives |G| < 1 always but this form has not been tried in the relaxation program. 
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C PPGGJ^AM TO CALCULATE THE A MPL I F I C A T I CN FACTOR GIVEN BY 

C FORMULA (E12) USING ALPHA GIVEN BY (E14). 

C 

C 

COMPLEX LI «ENU ,EMU , EMNU <• EMMU »T I «GNUy»GDEN 

RFAL»4 MU«NU 

IKN = 0 

dmax=o . 

R= 1 . 

R2=R*R 

N=5 

M = N + l 
DMACH=5./N1 

AL = 1 . 0 1 

Z I = ( 0 . , 1 . ) 

6 XM=l.00l 

_5 XM2- XM*XM 

02=7./{l.+5./XN?> 

DQ2=Q2/N 
!<N = 0 

KKN=0 

CMAX = 0 . 

C=SGRT ( C2) 

U2--D02 
CNEM= I .-XM.? 

■4 U2=U2 + 0Q2 

V2=Q2-U2 
U=SQRT ( AES ( U2 ) > 

V=SORT ( AES ( V 2 ) ) 

LiV=U*V 
D N U= P I / N 
D M U= P I / N 

MU-1 ,E-4 

DO 2 J-l,Nl 
NU= 1 .E-4 
DO 1 T= I .N 1 

KKN = KKN-f 1 

AL 1=Q2*XM2 
ENU-CEXP (Z I*NU ) 

EMNU=1 ./FNU 

EMU-CEXP(Z1*MU) 

EMMU=l ./EMU 

T 1 = V2^ ( EMNU’i‘E.MNU-2 . ’♦^EMNU + l . ) +2. *UV*R* ( I ,-EMNU-EMMU +E MNU ^^EMMU ) 
X +U2*R2*( 1 .-EEMMU>.''EMMU-2.*HMMU) 

Tl=Tl «CNEM+U2* <ENU-l . )- 0 « 5 *U V *R*ENU =»' ( E M U-EMMU ) 

^-KruF^T r- a:l T?Tcr“5p«TF.i7MG - r ; i - st ewNU-rn i 

Tl=U2* ( EMNU- I . 1 + V2*R2* ( EMU-2 .+EMMU )-0 . 5*UV*R*E.MNU* C EMMU-EMU ) 
GDEN=-T1+ALI *( U «R ( I .-EMMU I +V >M1.-EMNU)) 

0 1-C AB5 ( GNUM ) 

D2=C ABS( GDEN > 

D=Dl/02-l . 

IFCD.GT.O.) IKN=IKN+1 
IF ( D .GT .DMAX ) CMAX=D 
Too FORM AT ( 1 lEl 2 .4 ) 

1 NU-NU+ONU 

2 MU=MiU + DMU 

IF ( U2.GT .Q 2- 1 .E-3) GO TO 7 

GO TO 4 

7 PRINT 103.XM.IKN.KKN.DMAX 

103 FORMAT!/* MACH Na«.F6.3.' NO OF G GREATER THAN 1".I12. 

X • TOTAL COMPUTATIONS OF G«.I12. 

X * MAX G-l • «E^12.4) 

XM-=XM + DMACH 

IF (XM.LT .1 0 .01 ) GO TO 5 

STOP 

ENC 






